h1a_lm = lm_robust(attitudinal_outcome ~ h1_treat, data = survey_df) %>% tidy()
h1b_lm = lm_robust(behavioral_outcome ~ h1_treat, data = survey_df) %>% tidy()         

h2a_lm = lm_robust(attitudinal_outcome ~ h2_treat, data = survey_df) %>% tidy()
h2b_lm = lm_robust(behavioral_outcome ~ h2_treat, data = survey_df) %>% tidy()

h3a_lm = lm_robust(attitudinal_outcome ~ h3_treat, data = survey_df) %>% tidy()
h3b_lm = lm_robust(behavioral_outcome ~ h3_treat, data = survey_df) %>% tidy()


primary_mods = list(h1a_lm, h1b_lm, h2a_lm, h2b_lm, h3a_lm, h3b_lm)

primary_hypotheses = 
  bind_rows(primary_mods) %>%
  mutate(val = str_detect(term, "treat")) %>%
  filter(val == T) %>%
  mutate(outcome_lab = case_when(str_detect(outcome, "att") ~ "attitudinal compliance",
                                 str_detect(outcome, "beh") ~ "behavioral compliance")) %>%
  mutate(treatment_lab = case_when(str_detect(term, "1") ~ "H1: Any Appeal",
                                   str_detect(term, "2") ~ "H2: Personal Appeal",
                                   str_detect(term, "3") ~ "H3: Religious Appeal")) %>%
  mutate(treatment_lab = factor(treatment_lab, levels = c("H3: Religious Appeal", "H2: Personal Appeal", "H1: Any Appeal")))


ggplot(data = primary_hypotheses) + 
  geom_hline(yintercept=0, lty=2) + 
  geom_pointrange(aes(x=treatment_lab, y=estimate, ymin = (estimate - 1.96*std.error), ymax=(estimate + 1.96*std.error))) + 
  geom_pointrange(aes(x=treatment_lab, y=estimate, ymin = (estimate - 1.65*std.error), ymax=(estimate + 1.65*std.error)), fatten=2, size=1.1) + 
  labs(x='', y = '') + 
  ylab(expression("Average Treatment Effect (" ~ beta[1] ~ ")")) +
  facet_wrap(outcome_lab~.) + 
  coord_flip() + 
  theme_bw() +
  theme(legend.title=element_blank(), text=element_text(size=20), 
        axis.title.y=element_blank()) +
  theme(text=element_text(size=12, family="Times")) + 
  theme(panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(size = .2)) + 
  theme(strip.background =element_rect(fill="white"))

ggsave("./outputs/main_figure_1.pdf", width = 6, height = 3)


rm(list=setdiff(ls(), "survey_df"))
